/* 
   indexg - Index whole genome

   Copyright (C) 2007 

   This program is free software; you can redistribute it and/or modify
   it under the terms of the GNU General Public License as published by
   the Free Software Foundation; either version 2, or (at your option)
   any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU General Public License for more details.

   You should have received a copy of the GNU General Public License
   along with this program; if not, write to the Free Software Foundation,
   Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

#include <termios.h>
#include <grp.h>
#include <pwd.h>
*/

#include <stdio.h>
#include <sys/types.h>
#include <getopt.h>
#include "system.h"
#include "fastaIO.h"
#include "common.h"
#include "indexg_rep.h"

#define EXIT_FAILURE 1

static void usage (int status);

/* The name the program was run with, stripped of any leading path. */
char *program_name;

/* getopt_long return codes */
enum {DUMMY_CODE=129
};

/* Option flags and variables */

static char * ifile = NULL;		/* --input */
static char *oname = "stdout";		/* --output */
static UBYTE2 window = 12;		/* --window */
static BYTE2 step = 1;		/* --step */
static FILE * fhifile; 			/* input file handle */
static FILE *ofile;

static struct option const long_options[] =
  {
	{"output", required_argument, 0, 'i'},
	{"window", required_argument, 0, 'w'},
	{"step", required_argument, 0, 's'},
	{"help", no_argument, 0, 'h'},
	{"version", no_argument, 0, 'V'},
	{NULL, 0, NULL, 0}
  };


static int decode_switches (int argc, char **argv);

/* encode a short sequence (<=16nt) to integer code */
char *
decode ( BYTE4 code )
{
  char * seq;
  int i;
  int temp;
  seq = (char *)malloc(sizeof(char)*(window+1));
  seq[window] = '\0';
  for (i=1;i<=window;i++) {
	temp = code & 3;
	seq[window-i]=char_decode(temp);
	code >>= 2;
  }
  return seq;
}

int
main (int argc, char **argv)
{
  int m,n;
  BYTE4 i,j;
  UBYTE4 char_num;

  UBYTE4 encoded;
  UBYTE4 revencoded;
  UBYTE4 mask;

  char next_char;

  UBYTE4 bit_num;

  BYTE4 *big_table;
  FastaRecord * records;
  FastaRecord * cur_record;
  char * cur_seq;
  UBYTE4 cur_len;
  char * text;
  char * seq;

  char * seqname;

  BYTE4 NA = -1;

  program_name = argv[0];

  m = decode_switches (argc, argv);

  if (argc<=m) {
	usage(0);
  }

  /* do the work */
  bit_num = 1<<(window*2);
  mask = bit_num-1;

  seqname = (char*) malloc(sizeof(char)*10);
  
  /* Build big table */
  big_table = (BYTE4 *) malloc(sizeof(BYTE4)*bit_num);
  for (i=0;i<bit_num;i++) {
	big_table[i] = 0;
  }
	
  for (n=m;n<argc;n++) {
	ifile = strdup(argv[n]);
	fhifile = fopen(ifile,"r");
	if (!fhifile) {
	  fprintf (stderr, "cannot open `%s' for reading\n", ifile);
	  exit(1);
	} else {
	  fprintf (stderr, "process %s ...\n", ifile);
	}
 
	/* read sequence */
	fseek(fhifile,0,SEEK_END);
	char_num = ftell(fhifile);
	rewind(fhifile);
	text = (char*) malloc(sizeof(char)*(char_num+1));
	fread(text,sizeof(char),char_num,fhifile);
	fclose(fhifile);
	/* parse the whole text, produce a chain of fasta records */
	records = parse_text (text);
	free(text);
	
	cur_record = records;
	
	while (cur_record) {
	  encoded = 0;
	  revencoded = 0;
	  cur_seq = cur_record->fasta->seq;
	  cur_len = cur_record->fasta->length;

	  /* first window mer segment */
	  j = window;
	  while (j!=0 && *cur_seq) {
		encoded <<= 2;
		revencoded >>= 2;
		next_char = *cur_seq;
		cur_seq++;
		if (next_char=='N' || next_char=='n') {
		  j = window;
		  continue;
		} else {
		  encoded |= char_encode(next_char);
		  revencoded |= ((char_revencode(next_char))<<((window-1)*2));
		  j--;
		}
	  }
	  if (j != 0) break;			/* eof */
	  encoded &= mask;
	  big_table[encoded] ++;
	  big_table[revencoded] ++;
/* 	  printf ("%s\n",decode(encoded)); */

	  
	  while (*cur_seq) {
		j = step;
		do {
		  encoded <<= 2;
		  revencoded >>= 2;
		  next_char = *cur_seq;
		  cur_seq++;
		  if (next_char=='N' || next_char=='n') {
			j = window;
			continue;
		  } else {
			encoded |= char_encode(next_char);
			revencoded |= ((char_revencode(next_char))<<((window-1)*2));
			j--;
		  }
		} while (j!=0);
		if (j != 0) break;			/* eof */
		encoded &= mask;
		big_table[encoded] ++;
		big_table[revencoded] ++;
		/* printf ("%s\n",decode(encoded)); */
 	  }
	  cur_record =  cur_record->next;
	}
	free_records(records);
  }
  
  for (i=0;i<bit_num;i++) {
	seq = decode(i);
	fprintf (ofile,"%s\t%d\n",seq,big_table[i]);
	free(seq);
  }
  fclose(ofile);

  
  /* calculate occurence number for every position */
  for (n=m;n<argc;n++) {
	ifile = strdup(argv[n]);
	fhifile = fopen(ifile,"r");
	if (!fhifile) {
	  fprintf (stderr, "cannot open `%s' for reading\n", ifile);
	  exit(1);
	} else {
	  fprintf (stderr, "assigning value for %s ...\n", ifile);
	}

	/* open file to write */
	oname = malloc(sizeof(char)*(strlen(ifile)+4));
	oname = strcpy(oname,ifile);
	oname = strcat(oname,".rep");
	ofile = fopen(oname, "wb");
	if (!ofile) {
	  fprintf(stderr, "cannot open `%s' for writing\n", oname);
	  exit(1);
	} else {
	  fprintf (stderr, "write binary file %s ...\n", oname);
	}

	/* read sequence */
	fseek(fhifile,0,SEEK_END);
	char_num = ftell(fhifile);
	rewind(fhifile);
	text = (char*) malloc(sizeof(char)*(char_num+1));
	fread(text,sizeof(char),char_num,fhifile);
	fclose(fhifile);
	/* parse the whole text, produce a chain of fasta records */
	records = parse_text (text);
	free(text);
	
	cur_record = records;
	
	while (cur_record) {
	  encoded = 0;

	  cur_seq = cur_record->fasta->seq;
	  cur_len = cur_record->fasta->length;

	  

	  fwrite (cur_record->fasta->id,sizeof(char),10,ofile);
	  fwrite (&window,sizeof(UBYTE2),1,ofile);
	  fwrite (&step,sizeof(UBYTE2),1,ofile);
	  fwrite (&cur_len,sizeof(UBYTE4),1,ofile);

	  /* first window mer segment */
	  j = window;
	  while (j!=0 && *cur_seq) {
		encoded <<= 2;
		next_char = *cur_seq;
		cur_seq++;
		if (next_char=='N' || next_char=='n') {
		  /* rep number = -1 */
 		  for (i=j;i<=window;i++) {
			fwrite (&NA,sizeof(BYTE4),1,ofile);
		  }
		  j = window;
		  continue;
		} else {
		  encoded |= char_encode(next_char);
		  j--;
		}
	  }
	  if (j != 0) break;			/* eof */
	  encoded &= mask;
	  fwrite (big_table+encoded,sizeof(BYTE4),1,ofile);
	  
	  while (*cur_seq) {
		j = step;
		do {
		  encoded <<= 2;
		  revencoded >>= 2;
		  next_char = *cur_seq;
		  cur_seq++;
		  if (next_char=='N' || next_char=='n') {
			/* rep number = -1 */
			for (i=j;i<=window;i++) {
			  fwrite (&NA,sizeof(BYTE4),1,ofile);
			}
			j = window;
			continue;
		  } else {
			encoded |= char_encode(next_char);
			j--;
		  }
		} while (j!=0);
		if (j != 0) break;			/* eof */
		encoded &= mask;
		fwrite (big_table+encoded,sizeof(BYTE4),1,ofile);
	  }
	  cur_record =  cur_record->next;
	}
	free_records(records);
	fclose (ofile);
  }

  free(seqname);
  exit (0);
}
/* Set all the option flags according to the switches specified.
   Return the index of the first non-option argument.  */

static int
decode_switches (int argc, char **argv)
{
  int c;

  ofile = stdout;

  while ((c = getopt_long (argc, argv, 
						   "o:" /* output */
						   "w:"	/* window */
						   "s:"	/* step */
						   "h"	/* help */
						   "V",	/* version */
						   long_options, (int *) 0)) != EOF)
    {
      switch (c)
		{
		case 'o':		/* --output */
		  oname = (char *)strdup(optarg);
		  ofile = fopen(oname, "w");
		  if (!ofile)
			{
			  fprintf(stderr, "cannot open `%s' for writing\n", oname);
			  exit(1);
			}
		  break;
		case 'w':
		  window = strtod(optarg,NULL);
		  break;
		case 's':
		  step = strtod(optarg,NULL);
		  break;
		case 'V':
		  printf ("indexg_rep %s\n", VERSION);
		  exit (0);
		  
		case 'h':
		  usage (0);

		default:
		  usage (EXIT_FAILURE);
		}
    }
  
  return optind;
}


static void
usage (int status)
{
  printf (_("%s - \
Count n-mer for  whole genome\n"), program_name);
  printf (_("Usage: %s [OPTION] <FILEs>\n"), program_name);
  printf (_("\
Options:\n\
  -o, --output               output file for the number of every n-mer sequence, default: stdout\n\
  -w, --window               window size, default: 12\n\
  -s, --step                 step size, default: 1\n\
  -h, --help                 display this help and exit\n\
  -V, --version              output version information and exit\n\
\n\
e.g.\n\
  %s -w 12 -s 1 *.fa\n\
Note:\n\
The Output Binary File is divided in blocks, the format of every block is:\n\
Bytes\tContent\n\
1-10\tChromosome Name\n\
11-12\tWindow size\n\
13-14\tStep size\n\
15-18\tLength of Chromosome\n\
19-end of block\toccurence number for every Position till the end(-1 is for NA), 4 bytes each\n\
\n\
Author:\n\
  Tao Liu: taoliu@jimmy.harvard.edu\n\
"),program_name);
  exit (status);
}
